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Abstract 

We present numerical solutions of the two-dimensional Navier- 
Stokes equations by two methods; spectral and the novel Lattice 
Boltzmann Equation (LBE) scheme. Very good agreement is found 
for global quantities as well as energy spectra. The LBE scheme is, 
indeed, providing reasonably accurate solutions of the Navier-Stokes 
equations with an isothermal equation of state, in the nearly incom- 
pressible limit. Relaxation to a previously reported "sinh-Poisson" 
state is also observed for both runs. 
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1 Introduction 



In recent years 0, 0, [|, lattice gas models have been developed for a num- 
ber of fluid and fluid-like systems. Interest in the study of these methods 
derives from both theoretical and practical motivations. On the one hand, 
lattice gases provide a novel perspective on complex physical systems, dif- 
ferentiating the macroscopic physical effects that are observable from the 
simplified microscopic properties that ultimately are responsible for what 
is observed. From a computational perspective, lattice gases have shown 
promise as an alternative approach for solution of fluid-like partial differen- 
tial equations, especially on parallel computers where the relative indepen- 
dence of the lattice nodes can be exploited for computational efficiency Lat- 
tice gas models of the Cellular Automaton (CA) type have been developed 
for many systems, including hydrodynamics JIJ, magnetohydrodynamics[|]] 
(MHD), multi-phase flows ||, and flows through porous media|| |7]. Gener- 
ally speaking CA models are plagued with noise, so that very large spatial 
grids must be used. Schemes to improve the situation are made difficult 
by complexity of collisions, difficulty in eliminating spurious modes, lack of 
Galilean invariance, and other problems. Nevertheless, some appealing re- 
sults have been obtained using CA fluid models. Lattice Boltzmann (LBE) 
met hods [| represent an improvement in terms of noise and have produced 
promising results in computations. Many of the previous demonstrations of 
lattice gas have shown clearly that complex physical fluid phenomena can 
be reproduced by these methods. Examples include wave propagation in 
hydrodynamics and MHD, vortex "streets" in viscous flows around obsta- 
cles, immiscible fluid effects, and others. However, most and perhaps all 
of these demonstrations have been confined to a qualitative verification of 
the physics, and have stopped short of showing that the lattice methods in 
fact provide an alternative method for quantitatively accurate solution of 
the fluid equations. 

Since the introduction of LBE methods there has been recognition that 
the LBE framework provides opportunity to eliminate some, and perhaps 
all of the fundamental problems in the lattice gas approach. Two of these 
improvements are the use of the "single time relaxation approximation" 



(STRA) collisions [p], 11 1, and the ability to introduce corrections to the 



pressure that modify the equation of state and eliminate spurious modes] 11 



12]. In addition, the LBE approach permits greater flexibility in implemen- 



tations in terms of lattice structure and dynamical "rules". 
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Here we present a study of an improved LBE method in two dimensional 
(2D) hydrodynamics, demonstrating that for a simple nonlinear shear layer 
problem, the LBE method is accurate and effective. Related study of the 
performance of a 3D LBE scheme has been presented by Chen et alfll"3 



The problem we choose to address is a simple 2D periodic shear layer, per- 
turbed by the addition of a low level of random "noise". This basic flow 
problem is of broad relevance to flow applications in geophysics, aerody- 
namics and space physics, and has been studied in laboratory situations 
and through numerical simulations. Although three dimensional effects are 
absent in this treatment of the nonlinear shear instability problem, and the 
2D physical phenomena are well known, the simplicity and familiarity of 
this problem makes it a good starting place for accurate demonstration of 
the LBE method. 

The physical effects we are interested in reproducing are: spectral trans- 
fer, energy and enstrophy decay, relaxation to the long time "maximum 
entropy state" at times shorter than viscous decay time. Although our pri- 
mary purpose is to examine the incompressible behavior, the improved LBE 
scheme also is seen to accurately provide information about the "nearly in- 
compressible" features of the dynamics, including waves and nearly incom- 
pressible pressure and density fluctuations. 



2 Lattice Boltzmann Method 

We adopt a numerical scheme appropriate to 2D hydrodynamics that is 
based upon the Lattice Boltzmann Equation, giving rise to the abbreviation 
"LBE" method. This approach to solution of fluid equations is based upon 
the kinetic equations associated with Cellular Automaton (CA) models for 
fluids [1,2, 3]. In the CA formulation, the Boltzmann equation does not enter 
into numerical implementations, since it is constructed solely to demonstrate 
that certain averaged functions of the lattice dynamics approach solutions 
of the fluid equations. The LBE method arises from the suggestion |J that 
a direct solution of these equations would provide an alternative numerical 
approach to computation, conceptually midway between the Boolean CA 
dynamics and the continuum fluid equations. 

The several types of LBE models proposed thus far share with one an- 
other the advantage, relative to the underlying CA model, of significantly 
reduced noise. However, it has also been recognized that the LBE approach 
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allows for simplifications to improve numerical efficiency, as well as improve- 
ments to "cure" problems that arise in the underlying CA models. These 
refinements have substantially improved the prospect of useful LBE com- 
putations. The first major LBE refinement was the recognition that the 
"exact" LBE collision integral is unnecessarily complex and numerically 
inefficient . The first idea for streamlining the evaluation of the collision 
integral was|| to linearize the exact Boltzmann form. Evidently such a sim- 
plification preserves the tendency to approach the desired local equilibrium 
microscopic state, which is already known (from CA theory) to lead macro- 
scopically to hydrodynamics. The only cost is a certain amount of departure 
of the distribution from what it would be in the CA case. However, since 
the departures from equilibrium are generally assumed to be small, this is 
not expected to produce discrepancies in the physical results. Expanding 
on the idea that the details of the collision operator need not correspond 
to the Boltzmann approximation to the exact CA rules, two groups nearly 
simultaneously offered the suggestion fllO|, [U], [L2| that the exact collision 
operator can be, in effect, discarded, provided that one adopts a collision 
operator that leads, in a controllable fashion, to a desired local equilibrium 
state. By a "desired" equilibrium, we mean (1) one that depends only upon 
the local fluid variables, which themselves can be computed from the actual 
values of the local distribution at a point, (2) one that leads to the desired 
macroscopic equations (e.g., the Navier Stokes equation), and (3) one that 
admits whatever additional properties that are sought, such as simplicity 
or removal of nonphysical lattice effects. Recent work has shown that prop- 
erty (2) can be maintained rather easily, even when the collision operator 
departs significantly from the form taken in the Boltzmann treatment of the 
CA. In fact such departures are desirable from the point of view of several 
factors of type (3). 

Chen et alJIIJ offered the first suggestion that one could use the single 
time relaxation, or STRA, collision operator for a MHD LBE method. Sub- 



sequently, a similar method fTT|| was described, and referred to as a "BGK" 
collision integral, in reference to the more elaborate collision treatment of 
Bhatnagar, Gross and Krookflnj. The essence of the suggestion for the 



LBE method is that the collision term £l(f) be replaced by the well known 
classical single time relaxation approximation 

f _ feq 

n(/) = -^-. 

T 

An appropriately chosen equilibrium distribution is denoted by f eq , which 
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depends upon the local fluid variables, and a lattice relaxation time r that 



controls the rate of approach to this equilibrium. Later, Qian et al[ll] and 



Chen et al ||12|| described a STRA method for hydrodynamics that incorpo- 
rates the form (1), but which also includes a reservoir of "stopped" particles 
that enter into the equilibrium distribution to prevent the particle distri- 
bution from "cooling" in regions of higher fluid speed. The latter problem 
had plagued earlier CA implementations of fluid models by giving rise to 
a velocity dependent pressure. An improper equation of state of this kind 



introduces nonphysical compressive effects [T5|, [LB], including spurious os- 
cillations, and incorrect pressure profiles in channel flows. These effects are 
completely eliminated, to all orders in the Mach number, by these "pressure 



corrected" LBE schemes [11, 12]. In contrast, multispeed CA models only 
partially correct the equation of state, by moving the velocity dependence 
of the pressure to higher order. Still further improvements to the method 
were realized when the stopped particle reservoir was parameterized in such 



a way ||17|| that the sound speed could be controlled, enabling higher Mach 
number flows, and in principle, shocks, to be computed with the STRA-LBE 
scheme. 

In the subsequent sections, we present results obtained with an LBE 
scheme for 2D hydrodynamics, that incorporates a number of the above 
described features. We use a square lattice with eight moving particle di- 
rections plus stopped "particles" [II], |TS| . In CA terminology, this "9-bit 
model" refers to a lattice dynamical system in which particles stream from 
nodes on the lattice to the nearest neighbor nodes at fixed speeds, expe- 
riencing collisions at each node, which modify the particle state, and on 
average drive the particle distribution toward equilibrium. Nearest neigh- 
bor nodes relative to a node at x are located at the face-centers x + c^, 
for a = 1,2,3,4, with = (cos (a — l)7r/2, sin (a — l)7r/2), and the ver- 
tices of the square centered about x, i.e., x + c^ 1 , for a = 1,2,3,4, with 
cj/ = v^cos (a — l/2)7r/2, sin (a — l/2)n/2). To move to the appropriate 
node during the streaming step, a particle in state I a moves with velocity 
c T a while particles in the state II a move with velocity c ! J . In "lattice units", 
the lattice side can be taken to be 5x = 1 and the lattice streaming time 
5t — 1, so that type / particles have unit speed and type II particles have 
speed y/2. 

Turning to an LBE treatment of the dynamics, we denote the moving 
particle distribution function by for k = I or 77 and a = 1, 2, 3, 4, while 
the component of the particle distribution referring to the stopped particles 
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(which do not stream) is designated by /q. Adopting a single time collision 
operator, along with the foregoing streaming rules, we arrive at a kinetic 
equation 

/ Q fc (x + c*,T + 1) - / a fc (x,T) = - Ja Ja (1) 

T 

where k = I or //. 

Our LBE dynamical system is completed by choosing the equilibrium 
distribution |TT|], 



4 n 3 2i 
fo = gPtl-^ 2 ] 

f a = ||[l + 3^u + ^u) 2 -^ 2 ] (2) 

fa = |[1 + 3C^.U + ^.U) 2 -^ 2 ] 

where the mass density p and fluid velocity u are defined by 

P = fo + Efi ( 3 ) 

k,a 

and 

P U = Y, C afa- (4) 

k,a 

Several fundamental properties of this LBE scheme can be readily demon- 
strated, based upon the choice of equilibrium and the kinetic equation (1). 
It is straightforward to show that the pressure p = C^p, where C s = 1/VS 
is the sound speed. This is a "pressure corrected" LBE scheme with an ex- 
act isothermal equation of state. Next, considering moments of the kinetic 
equation, expanded according to a multiple scale Chapman-Enskog proce- 
dure, we find that the long wavelength low frequency behavior corresponds, 
at leading order, to a fluid equation for the velocity field. In addition, if one 
invokes a low Mach number ordering, which allows an approach to incom- 
pressible behavior, one obtains, in first approximation, the incompressible 
Navier Stokes equations, 

— + v • Vv = Vj9°° + z/V 2 v, (5) 

at po 

where p°° is the incompressible pressure, and po is the conserved initial 
mean density. Likewise, the small time dependent density variations obey 
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a continuity equation 

f + V-( P v) = 0. (6) 

At next order in the Chapman Enskog expansion, making use of the STRA 
collision operator with constant relaxation time scale r, we find that the 
viscosity is v — (2r — l)/6 (r > 0.5), which has the computationally de- 
sirable property of being independent of the density. Further discussion of 
the approach to incompressibility is given in section 6, and some additional 
remarks concerning the viscosity and the physical interpretation of r are 
given in the Appendix. 



3 Shear Layer Simulations: Spectral and Lat- 
tice Boltzmann 

The idealized shear layer consists of uniform velocity reversing sign in a 
very narrow region. That corresponds to a vorticity uj = (V x v) 2 differ- 
ent from zero only in the region of the sheared flow, and it is, in the ideal 
situation, a delta function. Therefore, we generated our initial conditions, 
in a simulation domain that is a square box of side 2n, with a spectral 
representation of delta functions at y = tt/2 and y = 3n/2 (with opposite 
signs) for the vorticity, truncated to include the appropriate Fourier ampli- 
tudes with wavevectors k = 1 through 8. This configuration is steady in 
the absence of viscosity and, although the simulations are viscous, we add a 
perturbation to trigger the nonlinear terms of the Navier Stokes equation. 
To this end the velocity Fourier modes with 1 < k < 60 where excited with 
random phases and with an energy spectrum of k~ 3 for high k, and peaked 
at k — 3. This "noise" was such that added about 10% of the energy al- 
ready present due to the idealized shear flow. In fact, a lower noise level 
would be adequate to trigger the nonlinear dynamics, but a larger level 
was used for reasons that will be further explored in Section 5. Thus, we 
have initially E k = 0.5, = 5.738644 (Enstrophy), P = 0.1866146E + 04 
(Palinstrophy) and Q = 0.2375568-E + 07 (Q-enstrophy). Proceeding from 
this initial condition, the subsequent dynamics gives rise to a familiar set 
of phenomena associated with the two dimensional shear layer, including 
vortex layer breakup, vortex rollup and coalescence of like-signed vortices. 
These will be described for the spectral and LBE runs in the next section. 

Operationally, the spectral run is a standard type, familiar in turbu- 
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lence computations, using an Orszag-Patterson implementation of a Fourier 
Galerkin scheme. We use a periodic box of side In and resolution 256 2 , 
which allows to use a Reynolds number R = 10, 000, and still resolve the 
dissipation wavenumber. This method solves the equation for the vorticity 
w = (Vx v) 2 , 

du/dt + u- Vuj = z/VV (7) 

using a fast transform evaluation of the nonlinear couplings, along with 
appropriate procedures for removal of aliasing errors. The Reynolds number 
for the longest wavelength is ~ 1/V. Time integration is a second order 
explicit scheme, using fixed time steps of At = 1/1024. The characteristic 
time scale for the motion of the large scale eddies is estimated from the 
energy as T SP = L/V2E, for a characteristic length scale L. In view of 
the slow decay of the energy, we estimate T$p using the initial value of 
the energy, and the unit length associated with the longest wavelength in 
the periodic box. Thus, in simulation units of time, the large scale eddy 
turnover time is Tgp ~ 1. The system is evolved up to simulation time t = 
119. The spectral simulation begins with the specified Fourier coefficients 
that generate the initial data, and the complete set of vorticity Fourier 
coefficients are stored at later times for subsequent comparison with the 
LBE results. 

For the LBE run, we obtain the initial fields from the relation v = 
Vip x z, where z is the unit vector in the z direction, taking the stream 
function ip from the solution to V 2 ip = —u, which is algebraically solved in 
Fourier space using the appropriate value of oj from the spectral run. The 
initial density is set to a constant. These fields are then used to initialize 
the distribution function / to its equilibrium value, for these specified fields, 
using Eq.(2). After this initialization procedure, the LBE system is evolved 
by subjecting it to the sequence of streaming and collisions alluded to above. 

For the LBE simulation we used a 512 2 box. For 2D hydrodynamics 
turbulence an estimation of the dissipation length scale Ld can be obtained 
with Ld/Lo ~ i? -1 / 2 , where Lq is the energy-containing length. For Lq = 
512 we get Ld ~ 5, that is, the viscous dissipation mechanisms are effectively 
occurring in a scale of the order of five cells. A smaller run of 256 2 size was 
carried out; however effects typical of lack of resolution of the dissipation 
length were observed for this computation. Having chosen the appropriate 
lattice size and (u 2 ) = 0.04 for the LBE system, the relaxation parameter 
r is then fixed to give the proper viscosity value to obtain R = 10, 000, using 
the expression v = (2r — l)/6 (in lattice units). In particular to achieve an 
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LBE Reynolds number, i.e R — UL/v — 10, 000, we use 



1 



(8) 



R = 0.04 x x 



2tt (2t - l)/6 



arriving as R = 10,000 when r = 0.500977848. Although the physical 
Reynolds numbers will be time dependent, scaling with the characteristic 
fluctuating fluid velocities, we expect that these LBE parameters produce 
Reynolds numbers in the two types of runs that are within 10% in value. 

To be able to compare LBE and spectral runs we also have to relate 
the time units, from lattice convection units (i.e., time needed to propagate 
microscopic information from cell to cell) to large scale eddy turnover time. 
That conversion is done in the following way. Using characteristic lengths 
and velocities for both spectral and LBE systems, we can find a relationship 
between the characteristic times for the schemes. Thus, using Llbe = 512, 



L SP = 2tt, Ulbe = y (u 2 ) = 0.04, Usp = 1, we can get an expression 
connecting the typical time for evolution of both systems, 



so we need about 2037 LBE time steps to complete one spectral character- 
istic time. 

Global features of the evolution are calculated dynamically for the LBE 
computation every 200 LBE-time steps (i.e. about 1/10 eddy turnover 
time). The velocity field is scaled to the units used for the spectral simula- 
tion and is Fourier transformed. The vorticity Cj{k) = i(k x v(k)) 2 is then 
evaluated, from which the energy, enstrophy, palinstrophy, q-enstrophy and 
mean square stream function are computed. This scheme assumes that 
the fluid is incompressible, i.e., for these global diagnostics, and for energy 
wavenumber spectra, the (small) admixture of nonvortical velocity fluctua- 
tions associated with the compressible LBE dynamics, is ignored. In Section 
6 we will further discuss the validity of this approximation. 




Tlbe _ Llbe/Lsp _ 512 1 
~%p~ ~ U LB e/Usp ~ ^TblM' 



(9) 



thus 



Tlbe = 2037.12 T SP 



(10) 
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4 Comparison of Spectral and LBE results 



In spite of the organized large scale appearance our initial data - a periodic 
shear layer perturbed by random fluctuations, the evolution of the system 
in time is quite typical of 2D incompressible hydrodynamics. Thus, the 
time histories of global quantities, illustrated in Fig. 1, are familiar in 
appearance and interpretation. Fig. 1 shows the evolution of the energy 
E = l w (k)| 2 /fc 2 , the enstrophy f2 = J2k ^(k)] 2 , the Palinstrophy P = 
J2k k 2 \u(]t)\ 2 , and (lacking a better nomenclature) the "Q-enstrophy" Q = 
J2kk 4 \u}(k)\ 2 . (The sums are over the independent wavevectors k.) 

E and Q are inviscid invariants, and therefore are monotonically de- 
creasing in these dissipative simulations. On the other hand P and Q can 
be amplified as well as dissipated, and are not monotonic. One can also 
prove that Q/E decreases in time ||1 9|| . This is associated with the tendency 
for 2D Navier Stokes flow to engage in "selective decay", wherein the tur- 
bulence drives the spectrum towards its geometrically determined extremal 
state fl = K^ in E, where K min is the lowest allowed value of wavenumber, 
in a time short compared with the decay of the flow due to viscosity. The 
perspective provided by Fig. 1 is consistent with prior results in showing 
that selective decay is at least a qualitatively useful picture of 2D turbu- 
lence. For example, focusing on Fig. la) and lb), one sees that E decays 
quite slowly compared with Q, allowing the conclusion to be drawn that the 
energy is "back-transferred" in k, where dissipation is slow. On the other 
hand, the flow tends to produce additional amounts of P (see Fig. lc), an 
effect that accelerates the dissipation of Q, since Q = —2vP. In addition 
Q is also amplified early in the run (Fig. Id), and is dissipated at later 
times along with E, Q and P. This complex process of spectral transfer, 
involving both direct transfer to higher k, and backtransfer to lower k is 



familiar in 2D flows p0[ , and is a consequence of a very large number of 
nonlinear couplings each involving triads of wave vectors. These couplings, 
as well as their symmetries that give rise to inviscid conservation of E and 
Q, are accurately simulated by the spectral method simulation technique. 
What is new in the panels of Fig. 1 is evidence that the LBE method tracks 
the spectral method closely with respect to evolution of E, Q, P and Q. 
Therefore, even though the LBE method does not involve a wavevector rep- 
resentation, or even the vorticity, in any direct way, it evidently provides 
a representation of the Navier Stokes dynamics that is accurate enough to 
preserve the subtleties of 2D nonlinear spectral transfer. 
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Each of the quantities E, f2, P and Q, provides a measure of the dis- 
tribution of vorticity over wavenumber, and those with higher powers of k 
weight the short wavelengths more heavily. A careful inspection of Fig. la-d 
shows that the quantities that emphasize the lower k part of the spectrum 
are most similar in the LBE and spectral runs. Evidently, the departures of 
the LBE from the incompressible spectral method are greatest at the higher 
wavenumbers. Nevertheless, even fine features of the spectral method evo- 
lution of P and Q are also seen in the LBE curves for the same quantities. 

Wavenumber spectra of the energy are compared in Fig. 2, for the spec- 
tral and LBE results, at times t — 0, 5, 49 and 80 (in simulation time units, 
i.e., eddy turnover times computed in terms of the initial data). Fig. 2a) 
shows, for the two runs, the initial spectra, which are identical by construc- 
tion. Local peaks at the lower wavenumbers are associated with the initial 
shear layers, while the higher k powerlaw is due to the "noise" perturbation. 
By t = 5, a substantial amount of spectral evolution has occurred in both 
LBE and spectral runs, but, as is shown in Fig. 2b), the energy spectra 
for the two cases have remained extremely close. The most noticeable de- 
partures are at the highest values of k, as expected from the discussion in 
the previous paragraphs. Very similar energy spectra are also seen at much 
later time, as is illustrated in Fig. 2c) and 2d) at times t — 49 and t — 80. 
In these latter two comparison plots one can see clearly that significant 
amounts of back transferred energy persists in the longer scales at these 
late times, and that this effect is accurately portrayed in the LBE run. 

Perhaps the most striking verification of the accuracy of the LBE run 
is found in the direct comparison of contour plots of the LBE vorticity 
with the spectral method vorticity at the same physical times. In Fig. 3 
we show pairs of vorticity contour plots at four times. While the times 
are given in simulation times, it should be noted that the equivalent LBE 
time was computed from the calibration discussed in the previous section. 
The early time state, at t — 1, is seen in Fig. 3a), which shows, in both the 
spectral and LBE cases, that the initial shear layers have begun the familiar 
process of vortex roll-up. The vorticity distribution is extremely similar in 
the two cases. By a later time (t = 5, in Fig. 3b) the roll-up has progressed 
and has produced individual vorticity concentrations. These subsequently 
convect in the flow due to all the vortices, and mergers occur between like- 
signed vortices due to "vortex collisions". Once again, the plots from the 
two methods show great similarity, even down to detailed structures near 
regions of like signed vortex interactions. Note that the same values of 
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vorticity contours are used in performing the comparisons. A distinctive 
vortex collision is captured by both methods at time t = 17, shown in Fig. 
3c. 

At later times, all the positive and negative vortex concentrations have 
separately merged into a single pair of vortices. Fig. 3d shows this state, 
computed in both the spectral and LBE runs, at t — 80. It is clear that the 
LBE method has succeeded in reproducing many of the important dynami- 
cal features obtained by the spectral method, which has been the standard 
method for turbulence. These features include the evolution of bulk quan- 
tities, the form and evolution of the wavenumber spectra, and the detailed 
features of vorticity contours, including vortex rollup and subsequent merg- 
ers of like signed vortices. We now turn to some more subtle features of the 
flow, which appear also to be well represented by the LBE method. 



5 Relaxation to "sinh-Poisson" most proba- 
ble state 

An interesting by-product of the decaying turbulence computation just de- 
scribed concerns the extent to which the two- vortex quasi-steady final state 
has vortex shapes which coincide with those recently seen at the end of a 
study of decaying two-dimensional turbulence reported elsewhere. A slight 
digression is required before it is possible to display the relaxation products 
of the turbulent computation in a way that will make this connection clear. 
It has long been realized that in decaying 2D Navier-Stokes flow, enstro- 
phy or mean-square vorticity decayed rapidly compared to energy or mean- 
square velocity, for reasons that are well known. The separation of the time 
scales increases with Reynolds number, and had led to a conjecture that the 
relaxed state of decaying 2D Navier Stokes turbulence would be one in which 



the enstrophy was minimized relative to the remaining energy [|T^, [21], |22|. ^3 
In such a state, the only excitations left in the energy spectrum would be 
those in the longest wavelengths allowed by the boundary conditions. Quali- 
tatively, such a "selectively decayed" state would resemble, for example, the 
states shown in the two panels of Fig. 3d. Some time ago|24], |25| , a highly- 
resolved (512x512), high Reynolds number (14,286, based on the largest 
allowed wavelength), and long-time (about 400 large-scale eddy turnover 
times) 2D Navier Stokes spectral- method computation was carried out, in 
an effort to test the above- described "selective decay" hypothesis. In broad 
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outline, the tests confirmed the hypothesis, but examined closely, departed 
from it. In particular, a scatter plot of computed pointwise vorticity vs. 
stream function revealed not a linear proportionality between the two, as the 
selective decay hypothesis would suggest, but rather a hyperbolic-sinusoidal 



one, in which the observed connection ||26|| was that the late-time vorticity 
and stream function were related by 

cuj = sinh(|/3|V0 (11) 
where c and (3 are constants. The result was surprising, to the extent that it 



had been predicted two decades ago [27, 28] from a mean- field theory of most 



probable states, not for a viscous Navier-Stokes continuum, but rather for 
a large number of ideal, discrete, parallel line vortices. The subject had de- 
veloped, with both analytical and numerical solutions of the "sinh-Poisson" 
partial differential equation |27|, |28| that had been derived to describe the 



most-probable, or maximum-entropy, states, and a good bibliography is 
given by Smith|2P|]. A reformulation of the maximum-entropy theory had 



been given in the context of magnetohydrodynamics[pO|, 3jJ , and a further 



development of the foundation of the Navier-Stokes basis for it will be given 
elsewhere |32|1 . Our intent in this Section is simply to point out that even 



this somewhat unexpected and perhaps exotic hyperbolic-sine connection 
between stream function and vorticity has been reproduced accurately in 
the present LBE computation. In Fig. 4, we graph two correlation functions 
vs. time, with the broken line referring to the spectral method computation 
and the solid line to the LBE computation. Shown in Fig. 4 are correlations 
between vorticity and stream function (lower curves) and between vorticity 
and the hyperbolic sine of j3 times the stream function (upper curves), where 
the constant j3 is determined from a least-squares fit to the computed data. 
For any two functions f(x,y) and g(x,y), the correlation C(f,g) is defined by 

r(fn)= (if - (f))(a - (g))) (u) 
u,9) l((f-(f)) 2 )((g-(g)) 2 )} 1/2 1 ] 

where the angle brackets () denote a spatial average over the entire box. 
Thus for any two functions which are proportional, C will be equal to unity. 
The approach to the "sinh-Poisson" prediction is seen not only to be far 
superior for the computed data, but it will also be noticed that the LBE 
and spectral method computations again track each other to a remarkable 
extent. We remark also that there is a problem, as yet unsolved, of extract- 
ing the observed statistical mechanical distribution of the LBE variables for 
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a vortex distribution directly from the LBE dynamics, without the neces- 
sity of detouring through the Navier-Stokes approximation. This must at 
present stand as a challenge for theory; a solution would be highly desirable 
as a logical link between the microscopic and macroscopic dynamics. 



6 Nearly incompressible hydrodynamics in 
the LBE scheme 

In the previous sections, evidence was presented that the LBE method repro- 
duces many of the essential dynamical features of the incompressible Navier 
Stokes equations, as computed by a spectral method code based upon the 
vorticity equation. In particular, we found that the solutions appear to be 
quantitatively similar to one another. What differences there are between 
the spectral and LBE results appear to be most pronounced at the higher 
wavenumbers. It is tempting to assign these discrepancies to "error" in the 
LBE formulation, and conclude that the methods correspond well, for most 
of the diagnostics of interest, at times of up to tens, or perhaps a hundred 
or more, characteristic nonlinear times. 

However, there remains the possibility that the LBE scheme, is, in some 
sense, more accurate than this would suggest. We refer here to the possi- 
bility that the LBE scheme, in effect is solving a compressible set of fluid 
equations, and therefore, would be expected to approximate solutions of the 
incompressible equations only in an appropriately defined limiting sense. In 
fact, the compressible Navier Stokes equation itself also possesses this prop- 
erty. For suitably chosen initial data, and for small Mach numbers, the 
solutions of the compressible fluid equations are expected |33|] to approxi- 
mate the solutions to the incompressible equations for at least some finite 
time interval. Simulations |34| of the compressible equations of 2D hydrody- 



namics (with a polytropic equation of state) have also led to the suggestion 
that finite Reynolds number extends the realm of this expectation, so that 
in some cases the "nearly incompressible" nature of a decaying flow may 
persist permanently. Since the LBE method is intrinsically compressible, 
we can reasonably expect that it, too, will admit a range of parameters and 
time in which its solutions approach the desired incompressible solutions. 
This, indeed, is what we have seen in the previous section. However, there 
is also the prospect that some of the departures of the LBE solutions from 
the spectral method incompressible solutions might be attributable to the 
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slight effects of compressibility. In that case, at least some fraction of the 
differences between the spectral and LBE solutions might not be errors of 
numerical origin, but rather physical effects that lie outside the realm of 
the incompressible equations. We briefly explore this possibility here, by 
examining whether the LBE results are consistent with the expectations of 
nearly incompressible fluid theory . 

Equation (5), the incompressible equation for the velocity field, assumes 
that the leading order velocity field, say, v , is divergenceless, V • v = 
and the leading order density is constant, say p = po = const.. Again 
assuming that this limit to incompressibility is obtained, we find that the 
incompressible pressure p°° appearing in (5), must satisfy 

VV° = -p V • (v • Vv ) (13) 

which is a consequence of the time independence of V • vo = 0. On the 
other hand, prior to the limit to incompressibility, the LBE system is found 
to obey the compressible equations 

p[^ + v- Vv] = -Vp + D, (14) 

which, along with the continuity equation (6) and the equation of state 
p = CgP completes the specification of the long wavelength, low frequency 
LBE dynamics. The term D on the right side of Eq. (14) represents the 
viscous dissipation terms for the compressible fluid limit of the LBE method. 
Since neither the form nor the effects of dissipative terms are central to the 
description of near-incompressibility that we examine here, we neglect D in 
the following discussion. 

Let us define the Mach number of the flow as M = 5v/C s , with 5v 
the rms value of v. When M << 1 we expect there to be conditions 
in which a decaying flow will remain nearly incompressible. Klainerman 
and Majda|33| have shown that the additional required conditions are that 



the initial data satisfy (|V • v] 2 ) 1 / 2 = O(M) and 5p = ((p - p ) 2 ) 1/2 = 
0(M 2 ) where (...) denotes a volume average. In the LBE run discussed 
above, the initial 5v = 0.04, and C s = l/v3, so the initial M = 0.069. In 
addition, p = p is uniform in the initial data. As for the velocity field, it 
is computed for the LBE initially in terms of the real space values obtained 
from the spectral method initial data. Thus, except for possibly errors due 
to the finite LBE grid, it satisfies V • v = initially. Consequently, the 
conditions for the approach of the compressible equations to the solutions 
of the incompressible equations appear to be well fulfilled. 
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In this circumstance, we expect that, for a finite time, the density should 
remain ordered as p = po + M 2 (p°° + p') + 0(M 3 ), while the pressure (in 
convection speed units) should satisfy p = M~ 2 (p + M 2 (p°° +p') + 0(M 3 )). 
Here, p°° is the incompressible pressure, satisfying the Poisson equation (13). 
There is also an additional pressure fluctuation p', at the same order as the 
incompressible pressure, but associated with acoustic waves, and decoupled 
from the incompressible equation of motion. The leading order density 
fluctuation is 5p ~ Af 2 (p°° + p') = 5p°° + 5p', where p' is also associated 
with acoustic waves. In addition to the Poisson equation, the incompressible 
pressure satisfies the relation p°°+p' « C 2 5p. In order for the incompressible 
dynamical equation to lack acoustic time scale variations, we must apportion 
the leading order density fluctuations so that p°° = C 2 8p°° = M~ 2 5p°°, the 
latter equality making use of the convection speed units. Considering also 
the velocity field, we note that, in a Fourier decomposition, we can readily 
divide the velocity field as v = v L + v_|_ where the longitudinal velocity v L 
has V • vl 7^ but V x = 0, while the transverse velocity vj_ satisfies 
V • Vi = but V x vj_ 7^ 0. Then, for maintaining near-incompressibility 
we require that the solutions remain ordered so that = O(M) for 
(incompressible) convection speed units in which vi = 0(1). 

The degree to which these expectations of nearly incompressible fluid 
theory are seen in the LBE solution can be examined by analysis of the 
LBE velocity and density fields. The results are illustrated in the panels 
of Fig. 5. The LBE velocity field is Fourier transformed and decomposed 
into transverse and longitudinal components by projections relative to wave 
vector k. Then, the rms transverse velocity U±_ = J (| vj_ | 2 ) and the rms 

longitudinal velocity Ul = \J (|vl| 2 ) are computed. The relative magnitudes 
of U± and Ul during the simulation are shown in Figs. 5a and 5b. In 
Fig. 5a, the time history of U±/C s is shown. The value decreases slightly 
from its initial value, which is very close to the value of the Mach number 
M = 0.069 computed from the entire velocity field. Thus, we expect that 
the longitudinal velocity is small, and this is verified in Fig. 5b, which shows 
Ul/U± as a function of time. The latter ratio meanders about a value near 
0.010. Consequently, in convection speed units, it is clear that the condition 
Ul = O(M), required for the approach to incompressibility, is well satisfied. 

The density fluctuations may be decomposed as well, in accordance with 
nearly incompressible theory. First, we simply evaluate the total density 
fluctuation 5p, and compare its value to the mean density and the Mach 
number as time progresses. This is shown in Fig. 5c as the solid trace, 5p/ po 
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vs. time. We see that the magnitude of the root mean square total den- 
sity fluctuation is comparable to the expected value, of order M 2 ps 0.0048. 
Next, we decompose the density into the part associated with the under- 
lying incompressible flow, and the part associated with acoustic activity. 
Using only the transverse velocity field, we numerically solve Eq. (13) for 
the incompressible pressure p°°. The incompressible density fluctuation is 
computed as 5p°° = M 2 p°°. The root mean square value of 5p°° is plotted 
also in Fig. 5c, normalized to the mean density. Again the result is clearly 
0(M 2 ). Finally, we compute the density fluctuation associated with leading 
order acoustic effects through Sp' = 5p — 5p°° at each point in space. Com- 
puting the root mean square Sp' provides a measure of the degree of acoustic 
activity. This is illustrated as well in Fig. 5c, showing that this component 
of the density fluctuation also remains of 0(M 2 ), again in accordance with 
the expectations of nearly incompressible theory. 



7 Discussion and Conclusions 



In the above sections we have presented a detailed comparison of solutions 
to the two dimensional Navier Stokes equations obtained from a Lattice 
Boltzmann method and from a more traditional spectral method. The flow 
problem considered was a familiar shear layer initial value problem, in pe- 
riodic boundaries and prepared initially with a low level of random noise. 
We find that the LBE method has provided a solution that is "accurate" in 
the sense that time histories of global quantities, wavenumber spectra, and 
vorticity contour plots, are very closely similar to those obtained from the 
spectral method. While the comparison is best at early times, the solutions 
remain extremely close to one another for at least several eddy turnover 
times, and in some ways remain close for times up to a hundred turnover 
times. In particular, details of the wavenumber spectra at high wavenum- 
bers are reproduced, as well as the detailed structure of vortex distributions 
seen in the contour plots. In addition, the LBE scheme faithfully repro- 
duces the recently reported long time tendency for the stream function to 
approach a "sinh-Poisson" state that emerges from a maximum entropy 
argument. We have also explored the possibility that the LBE solution, 
to the extent that it departs from a pure solution of the incompressible 
equations, is remaining in the mathematically delineated regime of "nearly 
incompressible flow" . This indeed appears to be the case, although a more 
complete verification would require comparison with a fully compressible 
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spectral algorithm, a refinement we have not as yet undertaken. 

It remains to discuss the accuracy of the LBE scheme in a quantitative 
way. To do so we have computed several kinds of normalized differences be- 
tween the results of the two runs, which are interpreted (for the most part) 
as errors in the LBE method. The normalized errors in the bulk quantities, 
energy, mean square stream function and enstrophy, are computed, for ex- 
ample, as \E SP — E LBE \/E SP , and shown in Table I. (The suffixes SP and 
LBE refer to <fi computed from the spectral or LBE schemes, respectively.) 
The normalized total rms error, defined for the spatially dependent variable 
(j) as 

^.p^sB)* (15) 

This rms normalized error has been computed as a function of time for <fi 
taken as u>, vj_ or ip. In addition we have computed the kurtosis K((f>) = 
(4> 4 ) j (4> 2 ) 2 for 4> taken as iv, v± or ip. Error in the kurtosis is conveniently 
expressed as AK/K = \K$p — Klbe\/ Ksp, where the suffixes have the 
same meaning as above. In table I we give the values of these normalized 
errors at spectral method times t — 1, 10, 50 and 100. 







TIME 


Error 




1 


10 


50 


100 






0.00742 


0.04867 


0.35675 


1.12567 


e 


V 


0.02136 


0.13685 


0.38283 


1.21901 




UJ 


0.12751 


0.53799 


0.65278 


1.37693 




i> 


0.00097 


0.00623 


0.00984 


0.00161 


AK/K 


V 


0.00297 


0.01966 


0.03168 


0.08017 




UJ 


0.00237 


0.01245 


0.05960 


0.05869 






0.00043 


0.00957 


0.01676 


0.01843 




E 


0.00081 


0.00075 


0.00045 


0.00035 




n 


0.01252 


0.01053 


0.01500 


0.00689 



It is immediately apparent that, at any fixed time, and for most cat- 
egories of error analysis, the error in ip is smallest, and the error in u is 
largest. In keeping with our previous discussion of the comparison of the 
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spectra, this is associated with the fact that the fractional error in the higher 
wavenumber excitations are greater than that of the lower wavenumbers. 

It is also apparent that the errors in the kurtosis are much smaller than 
the total rms errors, for a given field. In addition, the error in the bulk 
energy is less than the error in the kurtosis of the velocity. In fact, the 
kurtosis errors remain small compared to the total rms error, especially for 
the vorticity. The reasons for this appear to be that the structures, and 
the distribution of structures in the LBE run remain quite close to their 
spectral method counterparts. However, the exact positions of the struc- 
tures become different in the LBE case, relative to the spectral case. This 
disparity appears first in the high wavenumber structures, and later on in 
the large scale structures, so that by t = 80 (see Fig. 3d) the large vortices 
that remain are not at the same locations in the two runs. Nevertheless 
the spectra remain very close (see Fig. 2). As with the spectra, the kur- 
tosis calculation is not sensitive to the position of structures, but only to 
their magnitude and shape, and, in a statistical sense, to the distribution 
of shapes. Evidently, the distribution of excitations, in both wavenumber 
and real space, remains relatively close for the two methods. The largest 
error appears to be in the position of the vorticity structures, and the large 
increase in the error at later times is associated with the progressive drift 
in position of the LBE relative to the spectral results. 

The origin of the drift in vortex positions, while bulk quantities, shapes 
and spectra remain fairly accurate, can be attributed to several possible 
causes. First of all, to compare the methods, we needed to reconcile the 
LBE timescale with the spectral (fluid) timescale. This was accomplished 
(see Section 3) in the present study by computing a conversion factor at 
t — giving the ratio of the characteristic time units, involving the rms 
fluid velocity fluctuation. The latter quantity changes in time, but this 
change would not produce a difference in the results of the two methods if 
the fluid kinetic energy and the enstrophy remained exactly equal for the two 
cases. However, there is a small difference in the energies and enstrophies 
(see Fig. 1 and Table I), and this causes a slight inaccuracy in the times 
at which we compared the results. As these "clocks" drift apart, so do the 
positions of the vortices at the times at which we compare them. This part 
of the positional drift may be operational in our study, rather than intrinsic 
to the differences in the numerical methods, and could, in principle, be 
reduced by a more sophisticated, and more difficult, analysis of the data. 
A second cause of the positional drift, is closely related to the first, but is 
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of physical origin. Specifically, we have argued in Sec. 6, that some of the 
small departures from incompressible flow in the LBE method may be a real 
physical effect, that of nearly incompressible flow, which the LBE represents 
reasonably well, but which is absent in the purely incompressible spectral 
run. The effects of the small amount of compressible flow may include 
differences in the decay of energy in the two cases, as well as differences in 
the position of vortices, even at the same physical time into each run. In this 
perspective, the positional drift, as well as other differences in the results 
of the two methods, may be attributable to compressive effects, and not 
numerical error. We note that both of these possible sources of differences 
in the methods, are expected to have greater influence on the total rms 
error than on the bulk errors or the kurtosis error. This is because each of 
them induce small changes in the effective times of a comparison. During 
this small time increment the position of vortices vary rapidly compared 
to changes in the spectra, or compared to changes in their shapes (except 
possibly at times of vortex collisions). The total rms error is extremely 
sensitive to exact positions of all structures in the simulation domain, even 
if the structures are otherwise accurately represented. 

We are led to the conclusion that the LBE scheme has matured to the 
point that it offers an alternative method for solving incompressible flow 
problems with reasonably high accuracy. In particular, the above error 
analysis suggests that the LBE approach gives relatively good results for 
bulk quantities such as energy, for wavenumber spectra and for measures 
of distributions such as kurtosis. Although contour plots show great sim- 
ilarity in spectral and LBE cases, there is, evidently, a growing drift in 
relative positions of vortex structures in the two cases. However, for tur- 
bulence calculations, the importance of exact positions of the vortices is 
rarely considered central, while spectra, energy decay rates, and statistics 
such as kurtosis are of great interest. Moreover, we find some indication 
that the scheme also offers quantitative information concerning the small 
effects of compressibility, including "pseudosound" density fluctuations as- 
sociated with the incompressible flow, and accompanying acoustic waves. 
As far as efficiency is concerned, we note that, for these resolutions and at 
the Mach number used, the 256 2 spectral run "costs" about 6 cpu minutes 
per characteristic time, whereas the 512 2 LBE run "costs" about 8 min 
per characteristic time on the San Diego Cray YMP. Thus, the LBE is of 
comparable efficiency, and may fare better than the incompressible spectral 
code in a parallel implementation. However, one should also note that the 
timings of a compressible spectral code would be expected to be about a fac- 
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tor of M 1 longer to resolve acoustic frequencies. Consequently, if the small 
compressible effects are required, the LBE may already be more efficient. 

The particular LBE model we have used is the product of several re- 
finements to the method. These include corrections to the pressure that 
enforce a particular (isothermal) equation of state, and the use of a sin- 
gle time relaxation procedure for handling the collisional approach to local 
equilibrium. Further refinements and extensions are also feasible as well. 
In particular, the pressure can, in principle, be further modified to include 
an independent temperature variable, so that a full ideal gas equation of 



state can be implemented. In addition, the method can be modified |L7 
to allow for higher Mach number flows, and even transonic flows, to be 
computed. However, this has not been attempted here, in view of our goal 
of comparison with an incompressible solution, approached through a low 
Mach number flow. 
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Appendix 

In developing the LBE theory it is of interest to understand the relationship 
the theory has to kinetic theory of ordinary gases, in addition to evaluating 
the computational method itself. In this respect the LBE method described 
herein possesses some properties that are unusual from the ideal gas kinetic 
theory perspective. Specifically, the present model is developed to arrive at a 
useful computational representation of incompressible flow, evidenced by the 
emergence of Eq. (5) at lowest order in the Chapman Enskog expansion, and 
also at leading order in a Mach number expansion. However, particularly in 
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view of recent efforts to employ related LBE methods to flows that may 
be strongly compressible, it is important to examine features of the method, 
such as the viscosity, when compressible effects are included. Although a 
complete examination of these effects has yet to be completed, we have 
noted the following disparity between the simple STRA LBE method and 
ordinary gas kinetic theory. 

In the kinetic theory of simple gases, the kinematic viscosity is expected 
to be dependent upon density, approximately as v = \ij p, where the molecu- 



lar viscosity p is approximately independent of density |36|, |37|. This scaling 
emerges because one finds that p oc pv t h^, where v t h is a thermal speed 
(roughly analogous to the LBE lattice streaming speed) and A is a colli- 
sional mean free path, related to a collision time r c by A = v t bJc- I n spite of 
what appears as an explicit linear dependence of p upon p, it is a familiar 
result that the molecular viscosity is nearly density independent because A 
(or, equivalently r c ) scales as oc 1/p. More precisely, on the basis of kinetic 
theory, molecular viscosity is independent of density for a fixed tempera- 
ture, a fact originally noted by Maxwell, and born out in standard kinetic 



theory calculations (e.g.,|3^|). However, when such calculations are carried 



out with a single time relaxation approximation to the collision operator 
(with relaxation time r c ), the correct scaling is obtained only by associating 
with t c an inverse proportionality with density. 

The STRA LBE method used here and elsewhere [[TO, |TT| , |T2"| , differs 
from the ordinary gas kinetic theory result in that the relaxation time has 
typically been chosen as a density independent constant. Consequently, 
there are features of the LBE viscosity that differ from the ordinary gas 
situation. Most importantly, the molecular viscosity p is not independent 
of density, essentially because the combination pr still depends on density. 
The molecular viscosity cannot be immediately "pulled through" spatial 
derivatives, divided by p, and renamed as the kinematic viscosity v. Instead 
there are also new terms that appear, all of which involve Vp. This changes 
the form of the compressible dissipation terms (D in Eq. 14) to something 
other than the precise form expected for a compressible ideal gas. However, 
these additional terms, according to the nearly incompressible flow theory 
reviewed in Sec. 6, involve two more factors of Mach number than do the 
"usual" terms in the viscosity. Thus, the added effects do not directly or 
significantly influence the incompressible flow component of the LBE in the 
nearly incompressible regime. 

These differences reflect the fact that in LBE theory, in contrast to ordi- 
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nary gases (as well as cellular automata QTJ, 0, |!|), the collisional mean free 
part is not determined by actual collisions that occur in the dynamics. In- 
stead the "collision rate" is determined by the selected relaxation parameter 
that controls the rate of approach to local equilibrium. This parameter r is 
externally controlled and is arbitrary within the bounds set by stability con- 
ditions for the LBE dynamical equation. Accordingly, the constant STRA 
collision operator is adequate, and perhaps also an efficient way, to compute 
incompressible or nearly incompressible flow with an LBE scheme. However, 
an improvement may be desirable for LBE schemes that are designed for 
higher Mach number flows that admit more effects of compressibility |35[| . 
In particular, the STRA model can be modified by choosing r = r p /p 
with To a constant time scale, po the mean density, and p the local value of 
density. This modification is expected to bring a compressible LBE scheme 
into closer agreement with the kinetic physics of an ideal gas, particularly 
with regard to the structure and value of the viscous transport coefficients. 
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8 Figure Captions 

Fig. 1 Time history of a) energy, b) enstrophy, c) palinstrophy, and 
d) the next higher order moment, q-enstrophy (~ k A u(k) 2 ). Con- 
tinuous line corresponds to the LBE simulation. Departures are 
noticeable for the higher moments only. 
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Fig. 2 Wavenumber energy spectra for times a) 0, b) 5, c) 49 and 
d) 80, for both spectral and LBE simulations. Continuous line 
corresponds to the LBE simulation. The spectra for t — are 
identical for both runs by construction. 

Fig. 3 Isovorticity contour plots for times a) 1, b) 5, c) 17, and d) 
80. Dashed lines correspond to negative values of vorticity. The 
values for the contours are the same for all cases. Strikingly similar 
features can be found for the LBE simulation as compared with 
the spectral simulation. 

Fig. 4 Correlation between u and ip, and between u and sinh(|/3|^>) 
as a function of time. Continuous line corresponds to the LBE 
simulation. 

Fig. 5. Near incompressibility of the LBE run. a) time history of the 
rms transverse velocity normalized by the sound speed U±/C s . 
This quantity remains approximately constant, and equal to the 
initial Mach number M = 0.069. b) Ul/U± as a function of time, 
where Ul is the rms longitudinal velocity. This ratio is clearly 
bounded by M, as required for approaching incompressibility. c) 
Density fluctuations divided by p as a function of time for the 
LBE simulation, p, p°° and pf correspond to the total density, the 
"incompressible" density, and density fluctuations associated with 
acoustic waves, respectively (see text). All fluctuations are 0(M 2 ) 
(M 2 = 0.0048), consistent with nearly incompressible theory. 



9 Table Captions 

Table. 1. Normalized differences between the spectral run and the LBE 
run for various quantities for t — 1, 10, 50 and 100. e is the total 
rms error, whereas A$/$ = |$ 5 p — $lbe\/$sp- Large differences 
in e are due mainly to a drift in vortex positions. Differences 
are significantly reduced for the two lower sections of the table 
that show errors in quantities that are independent of the exact 
distribution of vorticity but are, instead, sensitive to the shape of 
the vortices. 



27 



